% RED NEURONAL
%-----------------------------------------------------------------------------------------------
% PARAMETROS DE LA RED
%-------------------------------------------------------------------------------------------------
%  version = 1;
%  hideLayerNeuronalNum = 10; 			% Numero de perceptrones de la capa oculta de la red
%  eta = 0.2; 								% tasa de aprendizaje de la red
%  beta = 1;
%  windowSize = 3;							% tamaño de la ventana de entradas de la red
%  countEpoch = 250;						% cantidad de epocas a procesar
%  activation_func = 1;					% Si esta en 1 cargo las se�ales normalizadas para la funcion
%  										% exponencial, sino cargo las de la tangente hiperbolica
%  										
%  momentum_enable = 0;					% flag para habilitar el calculo del momentum

alpha = 0.9;							% variable para el calculo del momentum
										
%eps = 0.00001;							% epsilon para la comparacion de 2 doubles.
eps = 0.000001;
successOnes = 0;						% Cantidad de salidas en uno correctos (aprendidos) de una epoca
Ones = 0;								% Cantidad total de unos de una epoca
b = 0.05;								% tasa de decremento de eta
a = 0.005;								% tasa de incremento de eta;
umbral = 0.3;							% tolerancia para considerar una salida como correcta
PreviousError = 999;					% Variable inicializada, utilizada para almacenar el error de
                                        % la epoca anterior
goodEpocs = 0;							% Variable que indica cada cuantas epocas adaptar eta
error= [];								% Vector que almacenara el error cuadratico medio de una ventana
error_epo = [];							% vector que almacenara los Promedios de los errores de cada epoca
etaVec = [];							% Vector que almacena los eta de cada epoca
successOnesVec= [];						% Vector que almacena las salidas en 1 aprendidas de cada epoca
epsEta = 0.000001;
%---------------------------------------------------------------------------------------------------
% CARGA DE VECTORES DE SEÑALES
%---------------------------------------------------------------------------------------------------
if( activation_func == 1 )
	load expsignals.mat;				
	signalsVec = expsignals;
else
	load tansignals.mat;
	signalsVec = tansignals;
endif
%------------------------------------------------------------------------------------------------
% CREACION DE MATRICES ALEATORIAS DE LA RED NEURONAL
%-----------------------------------------------------------------------------------------------
% random weigth vector of the first layer
% make sure that values of w in order of 1/sqrt(windowSize)
% w_kj is the connection between the input units with the first level hiden units
w_kj = rand(hideLayerNeuronalNum, windowSize);
w_kj = w_kj ./ 10;

% mixWeight makes vector have negative values.
for i = 1:rows(w_kj)
	w_kj([i], :) = mixWeight(w_kj([i], :));
endfor
  
% W is the connection between the first level hiden units with
% the output single unit
W = rand(1, hideLayerNeuronalNum);
W = W ./ 10;
W = mixWeight(W);
%--------------------------------------------------------------------------------------------
% AGREGADO DEL UMBRAL A LAS MATRICES DE LA RED NEURONAL
%--------------------------------------------------------------------------------------------
%add threshold to each row of weight
function ret = addThreshold(m, th)
	t = zeros(1, rows(m));
	t = th + t;
	ret = [m'; t];
	ret = ret';
endfunction

w_kj_with_threshold = addThreshold(w_kj, rand);
W_ji_with_threshold = addThreshold(W, rand);
%--------------------------------------------------------------------------------------------
% MATRICES AUXILIARES PARA IMPLEMENTAR MEJORAS EN EL ALGORITMO
%--------------------------------------------------------------------------------------------
W_ij_aux = W_ji_with_threshold;						
w_jk_aux = w_kj_with_threshold;						
%----------------------------------------------------------------------------------------------
% CREACION DEL VECTOR DE PATRONES CUYOS PATRONES POSEEN EL TAMAÑO DE VENTANA ESPECIFICADO PREVIAMENTE
%------------------------------------------------------------------------------------------------
windowVec = {};
l = length(signalsVec);
ceros = 0;
k = 1;
cantVecWithOutputInCero = 5000;				% Indica la cantidad de patrones cuyas salidas seran 0
for j = 1:l
	if( length( x = getWindow( j, signalsVec , windowSize)) != 0 )
		if( getWishedOutput( x, windowSize ) == 1 )
			windowVec{k} = x;
			k = k + 1;
		else	
			ceros = ceros + 1;
			% Si pasa el limite de ceros especificado no cargo mas vectores con salida en cero
			if( ceros < cantVecWithOutputInCero )
				windowVec{k} = x;
				k = k + 1;
			endif	
		endif
	endif
endfor
numPatterns = rows( windowVec' );					% cantidad de patrones
%-------------------------------------------------------------------------------------------------
% COMENZAMOS A CONTAR EL TIEMPO DE ENTRENAMIENTO
%------------------------------------------------------------------------------------------------
t1 = clock();
%------------------------------------------------------------------------------------------------
% ALGORITMO DE ENTRENAMIENTO
%-----------------------------------------------------------------------------------------------
for epoch = 1:countEpoch
	% valor del delta del peso del pattern anterior,
	% dif_W_ji_last y dif_w_kj_last son usados para el calculo del momentum
	dif_W_ji_last = 0;
	dif_w_kj_last = zeros(1, hideLayerNeuronalNum);

	t_ini_epoch = clock();

	for i = 1:numPatterns
		window = getWindowValues( windowVec{i}, windowSize );
		window = [window, -1];

		h_j = w_kj_with_threshold * window';
		
		h_j = h_j';

		v_j = zeros(1, hideLayerNeuronalNum+1);
		v_j = [g( h_j, beta, activation_func), -1];

		h_i = W_ji_with_threshold * v_j';
		output = g(h_i, beta, activation_func);

		wishedOutput =  getWishedOutput( windowVec{i}, windowSize );
		diferencia = abs( wishedOutput - output );
		
		if( wishedOutput == 1 )
		    Ones = Ones + 1;
			if( diferencia < umbral )
				successOnes = successOnes + 1;
			endif
		endif
		%Error cuadratico de este pattern
		err = 0.5 * ( (wishedOutput - output)^2 );

		%Guardo el error de este patron en un vector
		error = [ error, err ];

		deltaSubI = ( wishedOutput - output ) * gprima(beta, h_i, activation_func);
		if (momentum_enable == 1)
			dif_W_ji = (eta * deltaSubI * v_j) + alpha * dif_W_ji_last;
			dif_W_ji_last = dif_W_ji;
		else
			dif_W_ji = eta * deltaSubI * v_j;
		endif
		
		

		deltaSubJ = (gprima( beta,  h_j, activation_func) *  deltaSubI .* W);
		dif_w_kj = [];
		for k = 1:hideLayerNeuronalNum
			momentum = 0;
			if (momentum_enable == 1)
				momentum = alpha * dif_w_kj_last(k);
			endif
			dif_w_kj = [dif_w_kj; eta * window .* deltaSubJ(k) + momentum];
		endfor
		dif_w_kj_last = dif_w_kj;

		W_ji_with_threshold = W_ji_with_threshold + dif_W_ji;
		W = W_ji_with_threshold([1:rows( W_ji_with_threshold )], [1:columns(W_ji_with_threshold)-1]);
		
		w_kj_with_threshold = w_kj_with_threshold + dif_w_kj;
		w_kj =	w_kj_with_threshold( [1:rows( w_kj_with_threshold )], [1:columns(w_kj_with_threshold)-1]);
		
	endfor
	%-------------------------------------------------------------------------------
	% MEZCLA DE LOS PATRONES
	%-----------------------------------------------------------------------------
	windowVec = mixCell(windowVec);
	%--------------------------------------------------------------------------------
	% Mejora al algoritmo, adaptando la tasa de aprendizaje
	%-----------------------------------------------------------------------------------
	% Promedio del error de la epoca, formado por los errores cuadraticos medios de cada patron		
	etaVec = [etaVec, eta ];
	avgError = sum(error)/length(error);

	if (version == 1)
	% version == 1 es la version con las mejoras, en este caso se actualiza el eta.
		if(( PreviousError - avgError) < 0 )
			eta = eta - b * eta;
			W_ji_with_threshold = W_ij_aux;
			w_kj_with_threshold = w_jk_aux ;
		elseif(( PreviousError - avgError ) > 0 )
			PreviousError = avgError;
			goodEpocs = goodEpocs + 1;
			if( goodEpocs == 4 )
				eta = eta + a;
				goodEpocs = 0;
			endif
			W_ij_aux = W_ji_with_threshold;
			w_jk_aux = w_kj_with_threshold;
		endif
	endif

	%---------------------------------------------------------------------------------------
	error_epo = [error_epo, avgError];
	successOnesVec = [ successOnesVec, successOnes ];	
	printf("Error de la epoca %d: %f | tiempor transcurrido en la epoca %f [min]\n", epoch, sum(error)/length(error), etime(clock(), t_ini_epoch)/60);

	% Si la cantidad de salidas en 1 correctas es igual a las totales termino de entrenar
	if  ( successOnes == Ones )
		break;
	endif

	% verifico los ultimos 10 valores del error en cada epoca, si son iguales,
	% considero que se quedo estancado en un minimo local, entonces salgo del
	% programa
	if (length(error_epo) > 5)
		equals = 0;
		lastElement = error_epo(length(error_epo));
		for i = 1:4
			if (abs(error_epo(length(error_epo) - i) - lastElement) < eps)
				equals++;
			endif
		endfor
		
		if (equals == 4)
			break;
		endif
	endif

	if (eta < epsEta)
		break;
	endif

	%---------------------------------------------------------------------------------------
	% GRAFICOS  
	%-------------------------------------------------------------------------------------
	% Error de cada epoca en funcion de las epocas
	x = 1:1:columns(error_epo);
	subplot( 3, 1, 1);
	plot( x, error_epo );
	axis([1, length(x) + 1, 0, 0.05 ]);
	% tasa de aprendizaje ( eta ) en funcion de las epocas
	subplot( 3, 1, 2);
	plot( 1:length(etaVec), etaVec );
	axis([1, length(etaVec) + 1, 0, 0.3 ]);
	% Cantidad de salidas en 1 aprendidas en funcion de las epocas
	subplot( 3, 1, 3 );
	plot( 1:length(successOnesVec), successOnesVec );
	axis([1, length(successOnesVec) + 1, 0, Ones ]);
	drawnow;
	%-----------------------------------------------------------------------------------------
	printf("Cantidad de salidas en 1 totales = %d\n", Ones );
	printf("Cantidad maximo de 1 correctos = %d\n", max(successOnesVec) );
	printf("Cantidad de salidas en 1 correctos en la epoca actual = %d\n\n", successOnes );
	successOnes = 0;
	Ones = 0;
	error = [];
	fflush(stdout);
endfor

%---------------------------------------------------------------------------------------
% GRAFICOS
%-------------------------------------------------------------------------------------
print("grafico_entrenamiento.png");
%-----------------------------------------------------------------------------------------


%---------------------------------------------------------------------------------------------------
% RESULTADOS
%---------------------------------------------------------------------------------------------------
printf("\n------------------------------\nRESULTADOS\n");
if (version == 1)
	printf("Version mejorada\n");
else
	printf("Version normal\n");
endif

if (activation_func == 1)
	printf("Funcion de activacion: exponencial\n");
else
	printf("Funcion de activacion: tangente hiperbolica\n");
endif
	
printf("Cantidad maximo de 1 correctos = %d\n", max(successOnesVec) );
printf("Tiempo transcurrido: %f [min]\n", etime( clock() , t1 )/60 );
printf("Eta inicial: %f\tEta final: %f \n", etaVec(1), min(etaVec));
printf("Neuronas en la capa oculta: %d\n", hideLayerNeuronalNum);
printf("Tamanio de la ventana: %d\n", windowSize);


%---------------------------------------------------------------------------------------------------
% GUARDAMOS LAS MATRICES DE LA RED NEURONAL
%--------------------------------------------------------------------------------------------
save w_final.mat w_kj_with_threshold;
save W_final.mat W_ji_with_threshold;
%----------------------------------------------------------------------------------------------
